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ABSTRACT 

Young protostars embedded in circumstellar discs accrete from an angular momentum- 
rich mass reservoir. Without some braking mechanism, aU stars should be spinning at 
or near break-up velocity. In this paper, we perform simulations of the self-gravitational 
collapse of an isothermal cloud using the ORION adaptive mesh refinement code and 
investigate the role that gravitational torques might play in the spin-down of the dense 
central object. While magnetic effects likely dominate for low mass stars, high mass 
and Population III stars might be less well magnetised. We find that gravitational 
torques alone prevent the central object from spinning up to more than half of its 
breakup velocity, because higher rotation rates lead to bar-like deformations that 
enable efficient angular momentum transfer to the surrounding medium. We also find 
that the long-term spin evolution of the central object is dictated by the properties of 
the surrounding disc. In particular, spiral modes with azimuthal wavenumber m = 2 
couple more effectively to its spin than the lopsided m = 1 mode, which was found 
to inhibit spin evolution. We suggest that even in the absence of magnetic fields, 
gravitational torques may provide an upper limit on stellar spin, and that moderately 
massive circumstellar discs can cause long-term spin down. 

Key words: stars: rotation — stars: protostars — hydrodynamics — methods: nu- 
merical — accretion, accretion discs 



>. 1 INTRODUCTION 



One of the unsolved problems in the physics of s tar forma- 
tion is the spin of protostars (jBodenheimeil [l995l ) . The ob- 
served specific angular momentum of molecular clouds typ- 
ically excee ds that of stars by a s much as four orders of 
magnitude (jCoodman et al.lll993l ). Most of the modern at- 
tempts at explaining the removal of excess angular momen- 
tum from the prot ostar invoke magnetic t o rques or magnetic 
stella r winds (e.g. iMatt fc Pudrit j |2005| . l2008l : iMatt et all 
I2OI0I . and references therein). Although some theoretical 
studies on the effect of purely hydrodynami c star-disc inter- 
actio n s on stellar spin have been c a rried out (I Yuan fc Cassenl 
19851: IPopham fc Nara^jil Il99ll : iBisnovatvi-KoganI 1 199 J 
Glatzel fc ObachI 1 1999*) .~ non-magnetic mechanisms have 



largely been neglected. However, there may exist situations 
where magnetic fields are unavailable to remove angular mo- 
mentum. 

Population III stars have long been thought to be 
unmagnetised due to both the absence of seed fields, 
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and the high temperatures required to m aintain a sufR- 
cient ionisation fraction in metal free gas (|Tan fc McKed 
'2004L 

However, these ideas have recently b een challenged 
(Fe derrath et al.|[20^11 : Schleicher et al.ll201lh . and the issue 
remains unresolved. Recent simulations suggest that with- 
out the action of magnetic fields or gravitational torques, the 
first stars might reach near break-up velocities (|Stacv et al.l 

Even for present-day massive star formation it is un- 
clear that the magnetic mechanisms normally invoked to 
regulate stellar spins are applicable. Massive stars form 
with much higher accretion rates than low-mass ones, and 
their discs tend to be dominated by gravitational rather 
than magnetic angular momentum transport mechanisms 
dKm mholz et al. 2007; Krattcr et al. 2008; Krumholz et al] 
l2009l : lHennebelle et al.ll201ll : |Peters et al.ll201ll '). Naively in- 
serting the high accretion rates typical of massive star for- 
mation into the models mo st commonly adopted to explain 
T Tauri star spin rates (e.g. lMatt fc Pudrit j|2008l ') yields the 
con clusion that massi ve stars should be spinning at break- 
up (|Rosen et al.ll2011, in preparation) , in contrast with ob- 
served stellar spins I Wolff et al.ll2006l ). Recent observations 
of discs around young B stars are also consistent with non- 
magnetised accretion columns, in contrast to their lower 
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mass counterparts l|Eisner et al.ll201ol ). Thus, non-magnetic 
braking mechanisms are also of interest for hmiting stellar 
spin. Since the most general shape a protostar can have is 
triaxial, gravitational torques from a protostar acting on its 
surroundings should be explored. 

Gravitational torques as an important means of angular 
momentum transport have been proposed to solve the angu- 
lar mo mentum problem of cloud cores in an analytic work by 
iFisheJ 12004) and a numerical study bv ljappsen fc KlessenI 
1 2004 ). However, these authors consider multiple stellar sys- 
tems. Instead, the simpler problem we consider — that of a 
single object interacting wit h its surroundin g s — is closer 
to the setup envisioned by lYuan fc CassenI (|l985l V Using 
existing theories, they discussed how the spin angular mo- 
mentum of a protostar may be regulated by the circum- 
stellar disc. If the protostar rotates rapidly, it may be sub- 
ject to instability and become triaxial. It will then exert a 
torque on the disc, allowing the protostar to transfer angu- 
lar momentum to the disc. To th e best of our kno wledge, no 
numerical simulations similar to lYuan fc Cassenl 's problem 
have yet been undertaken. However, transport of angular 
momentum via gravitational torques has been observed in 
simul ations focusing on the first core phase of star forma- 
tion llBatelll998l:ISaigo fc Tomisakall2006l : ISaigo et al.ll2008l : 
ISaigo fc Tomisakall201li l 

In this work, we present calculations of the collapse 
of an isothermal sphere and study the role of gravitational 
torques acting on the central object to remove its spin angu- 
lar momentum. By modelling the central object as a finite- 
volume fluid body, we demonstrate through simple numeri- 
cal experiments that its spin can be limited by its deforma- 
tion and the gravitational interaction with the surrounding 
medium. In addition, we find that long term spin-down is 
possible if the surrounding disc develops m = 2 spiral modes, 
but spin-evolution is inhibited if the disc develops significant 
m = 1 non-axisymmetry. 

This paper is organised as follows. In 32]we describe our 
model setup and numerical method. We define diagnostic 
measures in iJS] We compare and contrast two simulations in 
331 pointing out important correlations between disc modes 
and spin-down. We discuss several important caveats to our 
conclusions in !j5l In !j6]we discuss implications of our results 
on protostellar spins and stellar evolution. 



2 STAR-DISC MODEL 

We consider the collapse of a cloud leading to the formation 
of a dense central object surrounded by a disc. For conve- 
nience we will use 'star' and 'stellar' to refer to the central 
object and its properties, even though the problem as we 
set it up is fully dimensionless, so the central object does 
not necessarily correspond to the physical size or internal 
structure of a real star. We discuss this issue further in i)2.2l 
The system has mass Msys — Md + M, , where Ma and 
M* are the dis c and stellar masses respectively. A recent nu- 
merical study jKratter et al.|[2"oiOl ) shows that discs formed 
via this idealised collapse are characterised by two dimen- 
sionless parameters describing its mass accretion and rota- 
tion rates. 

Accretion of material onto the disc from the cloud is 



described by the parameter ^: 



GM 



(1) 



where M is the mass accretion rate and Ciso is the isother- 
mal sound-speed. The cloud rotation, responsible for disc 
formation, is described by the parameter F: 



M 



(2) 



where flk is the Keplerian orbital frequency of material join- 
ing the system from the cloud, assumed to occur at cylindri- 



cal radius Rk such that 
that 



GMsys/ Rk- If can be shown 



(3) 



where Msys = Mt has been used, t is the elapsed time and 
h is the disc aspect-ratio at R^. In this work we specify h 
instead of F directly. 

The system is modelled as a single, non-magnetic, in- 
viscid and self-gravitating fluid with density p, pressure P, 
velocity v and gravitational potential <I>. Its evolution is gov- 
erned by the usual Euler equations and the Poisson equa- 
tion. To distinguish the central object or star, the pressure 
is calculated via a barotropic equation of state (EOS) 



2 71 
CisoP 



1 + 



(4) 



where 71 = 1.00001, 72 = 5/3 and p« is a fixed transition 
density defined below. 

This EOS mimics star formation by halting gravita- 
tional collapse. For p <gi pt, the fluid collapses isothermally. 
As the collapse proceeds, the cloud's central density in- 
creases. When p ^ pt, P x p"'^^ and further collapse is 
prevented by the increased pressure. The central object is 
then effectively a polytrope with polytropic index n = 3/2. 
Because it has finite volume, it can deform. The resulting 
non-axisymmetric object may then be spun down due to 
gravitational torques from the external, lower density disc. 

2.1 Initial conditions 

The system is initially an isothermal sphere of radius rc. 
For r ^ r« = qr^, where r is the spherical radius and q is a 
dimensionless parameter, the initial density is 



Po{r) 



Act 



(5) 



47rGr2 

The dimensionless parameter A relates to the a ccretion rat e 
^, and we use tabulated values of ^4-^ pairs from[Shul (|l97'i1 ). 
The spherical region r ^ r« is designated as the initial star. 
We set p* — po{r») and po{r < r,) = p«. 

The cloud is initialised with an azimuthal velocity 



W0 



R/r 



1 



P ^ r. 
R>r* 



(6) 



where R is the cylindrical radius. The initial star has solid 
body rotation and is below break-up speed at its equatorial 
plane but faster than the rest of the core. The latter may bias 
angular momentum loss from the star to the disc, because 
if the star becomes triaxial, the pattern speed of the stellar 
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potential felt by the disc may be naturally higher than the 
disc rotation frequency, exciting trailing spiral density waves 
in the disc. 

Apart from the small region r ^ qr^, the init i al con - 
ditions here are identical to that of iKratter et alj l|2010t ) . 
Our motivation for this setup is the same as in that paper: 
this setup produces a collapse in which the dimensionless 
parameters f and F remain fixed, and thus produce a par- 
ticularly clean numerical experiment for studying behaviour 
as a function of ^ and F. A real protostellar core, of course, 
will not follow this simple structure, and ^ and F will vary as 
the core accretes. However, this variation will generally be 
on time scale s that are lone compared to the disc or bital pe- 
riod (e.g. see lKratter et al'lboosf : lOffner et al.lboiol ') and we 
may therefore, to first order, think of the system as simply 
proceeding through a series of quasi-equilibrium states with 
fixed ^ and F. Thus experiments at given ^ and F provide 
valuable insight. 

It should be noted, though, that in the idealised collapse 
described by constant ^, F, the central star is a point mass 
that can neither gain nor lose angular momentum, whereas 
in our problem it can freely exchange angular momentum 
with the disc. This effect is neglected in simple collapse mod- 
els, so our cores may evolve away from it. With this in mind, 
it is best to regard 5, F as initialisation parameters. 



2.2 Dimensionality and physical scales 

As in the self-sim ilar collapse of IShul (|l977l ) and in 
iKratter et al.l l|2010t ). the problem we have specified is fully 
dimensionless, and so we report our results in terms of non- 
dimensional numbers. In the evolutionary plots shown be- 
low, time is scaled by r = 2Tvq{rc/ciso)\/S/A, which is the 
Keplerian orbital period at r = qvc- 

Because our problem is dimensionless, we can scale our 
simulations to apply to a range of physical systems, in the 
context of low and high mass star formation (and possi- 
bly even gaseous planet formation). However, there is one 
caveat: the dynamic range we are able to achieve in our sim- 
ulation is much less than the dynamic range involved in the 
formation of a real star (r ~ lO^^cm) from a real protostellar 
core (r ~ lO^^cm). Even with adaptive mesh refinement, we 
achieve a dynamic range of ~ 10*, not ~ 10^. If one wishes 
to scale our dimensionless experiments to physical scales, 
one may think of doing so in two ways that give dynamic 
range comparable to what we achieve. 

First, one could envision that we are studying the col- 
lapse of a protostellar core leading up to the formation of the 
first hydrostatic core that is ~ 5AU in size (|Masunae:a et al.] 
1 19981 ). Alternately, one could envision that our central object 
is the size of a star, but that we are simulating the evolution 
only of the inner ~ 10^^ cm of the core that surrounds it. 
However, there is no evidence that the physical effects we 
identify depend in the slightest on the dynamic range of the 
simulations. 



2.3 Numerics 

The hydrodynamic equation s are evolved using the 
Godunov-type OR ION code l|Truelove et al.1 Il998l : iKleinI 
ll999l : lFisheill2002l ') in Cartesian co-ordinates {x,y,z) . The 



computational domain is a cube of length L = Arc- ORION 
offers adaptive mesh refinement, a key advantage for the 
multi-scale flow considered here. We do not make use of 
ORION'S sink particle or radiative transfer capabilities in 
this study. We use a base grid of 128^ with 6 levels of refine- 
ment, giving the highest effective resolution of 8192"^. 

A characteristic length-scale for self-gravitating prob- 
lems is the Jeans length. 



(7) 



where Cs = \J dP/dp is the density-dependent sound-speed. 
For a grid spacing 5x, a measure of resolution is 5x/\,j. Our 
EOS gives 



2/3' 



1/2 



(8) 



For example, if p/p, ~ 10, then the linear resolution is a 
factor of ~ 3 better than it would be for gas of the same 
density obeying an isothermal EOS. 

We use a Jeans number Nj = 8 to define the maximum 
resolvable density pj 



PJ 



TV 

NjSx V G 



(9) 



and we refine \i p > p,j at each grid level. Note the above is 
the isothermal Jeans density, which is smaller than the true 
Jeans density based on our EOS. Thus, using pj as defined 
above is a conservative approach. To ensure the star-disc 
interface is resolved, we also refine if p > 0.5p*. We also 
refine to the highest level within the theoretical disc radius 
Rk, and within one scale- height in z. 



3 DIAGNOSTICS 

In this section we define the diagnostic tools used to inter- 
pret our results. Measurements of the star are summarised 
in Table [T] where integrals are taken over the region p ^ p« 
[dV is the volume element); symbols preceded by A are rel- 
ative to the star; 'KE' and 'PE' stand for kinetic and poten- 
tial energies, respectively. Note that <1> is the gravitational 
potential due to all of the fluid. 

We regard fluid with p ^ p, as stellar. The star has 
position and velocity (a;*,v«). The stellar spin angular mo- 
mentum is defined with position and velocities of fluid ele- 
ments with respect to (a;«,i;«). The stellar rotation radius 
5** is deflned from its moment of inertia and is used to de- 
fine stellar spin and break-up frequencies. 5", is not to be 
confused with the stellar surface where p = p, , which is typ- 
ically larger than 5** and generally non-axisymmetric. We 
use (r*) to denote the average radius in the star's equatorial 
plane where p = p* . 

The evolution of star's angular momenta is strongly in- 
fiuenced by the fluid external to the star. The properties of 
the fluid are measured in cylindrical co-ordinates [R, 4>, z) 
centred on the star with velocities relative to u«. We take 
the vertical direction of the cylindrical co-ordinates to be 
parallel to the vertical (z) direction of the inertial frame, 
assumed to be aligned with the stellar spin axis. The disc 
mass Md is defined to be the difference between the total 
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mass within a cylinder of radius i?^, thickness 2hRk, cen- 
tred about the star, and Af, . The disc mass is only defined 
when the mass within this cylinder exceeds M, . 

For simplicity, we integrate the external fluid vertically 
over a slab of constant thickness and work with surface den- 
sity S = J pdz and vertically- averaged velocity U, where 



Table 1. Definition of stellar measurements. 



u 



^1 p{v 



v»)dz. 



(10) 



Our results are insensitive to the extent of vertical integra- 
tion, and in fact also insensitive to whether we use relative 
velocities or those in the inertial frame. 

Non-axisymmetric modes in the surface density of the 
external fluid can be found via Fourier analysis defined by 



am{R) 



E(i?, (p) exp {~im- 



(11) 



where a^n (R) is the radius-dependent amplitude and m is the 
azimuthal wave-number. The integrated amplitude Cm = 
J UmdR is used as a global measure of mode amplitudes. 
This integration excludes the star. 

The relevant angular momentum fiuxes for this problem 



Fa = RT^U^Ur 
Fr = RT.5U4,5Ur, 

Fg^ [ dzdR^d^^/'inG. 



(12) 
(13) 

(14) 



Fa is the flux associated with large-scale advection and 
Fr is the flux due to Reynolds stresses, where S here de- 
notes deviations from azimuthally-averaged values. Fq is the 
angular momentum fiux associated w ith self-gravitational 
torques (|Lvnden-Bell fc Kalnaislll972l '). These fluxes can be 
converted into a viscosities denoted qa, ctR, olg e.g.. 



OLG 



Fa 



R 



dlnn 



din 7? 



(15) 



and similarly for a a , ctR ■ Q = U,/,/ R is the angular velocity 
and (•)</, denotes an azimuthal average. 

It is important to regard the a's above simply as non- 
dimensionalised fiuxes rat her than the viscosity coe fficient 
in accretion disc models (jShakura Sz SunvaevlllQT^ . which 
would impose ^ a < 1). A positive fiux out of the star 
indicates spin angular momentum loss. To ensure that this 
angular momentum loss is physical rather than numerical 
(since ORION does not exactly conserve angular momen- 
tum), we can compare the physical a's we measure to the 
effective numerical viscosity, denoted by Q iv, measured for 
the ORION code bv lKrumholz et all l|2004 ). 

It will also be useful to consider mass accretion rates. 
To measure the ^ parameter from our simulations, we define 



^sin 



-3--(M, + M. 



(16) 



Although this is the most natural definition, because fre- 
quent simulation output is not practical due to the large 
data sizes involved, the numerical time derivative only gives 
an estimate. Nevertheless, we find that the mass accretion 
rate is consistent with the behaviour of stellar spin. 



Name 


Symbol 


Definition 


Mass 




JpdV 


Position 




f xpdV/Mt 


Velocity 


Vt 


f vpdV/M, 


Rotation radius 




{/ [Ax2 -1- Ay2] pdV/M,Y'^ 


Orbital ang. mom. 


jo 


a;, A t), ■ z 


Spin ang. mom. 


is 


/ pAa; A Av ■ zdV/M* 


Spin freq. 






Break-up freq. 


Ob 




KE-to-PE 


T/\W\ 





4 RESULTS 

We study the evolution of the star's spin in two cases. In 
Case 1, we consider initial conditions leading to a large 
disc-to-star mass ratio. For our comparison (Case 2) run 
we consider initial conditions which lead to a smaller disc- 
to-star mass ratio. We shall see that on the overall simulated 
timescale, the massive disc (Case 1) is less efficient at spin- 
ning down the star than its lower mass counterpart (Case 
2) due to different azimuthal spiral modes dominating the 
angular momentum transport. 

Before discussing our two cases, we should men- 
tion an important constraint on our parameter choices. 
iKratter et ahl (|2010l ) have established the ranges in parame- 
ter space for which disc fragmentation is or is not expected. 
We limit ourselves to values of ^ and F such that the disc 
will be non-fragmenting within the simulation. We do so in 
order to render the numerical experiment as clean and easy 
to interpret as possible. 

Massive stars form with high accretion rates that 
do te nd to produce gravitationally-unstable, fragmenting 
discs llKrumholz et al.ll2007l.l2009l:lKratter fc Matzneij|2006l: 



Kratter et all l2008l: iPeters et al.1 l2010l: iHennebelle et all 



201ll ). However, protostellar discs are heated by the cen- 



tral star and by viscous dissipation, both of which increase 
at smaller radii, so that fragmentation generally occurs only 
outside ~ lOOAU. Since it is the angular momentum flux 
out of the star to the disc that counts, it is the inner disc 
that is important for star-disc angular momentum exchange. 
This means that our choice to limit ourselves to the non- 
fragmenting part of parameter space does not prevent us 
from applying our models to massive protostars and their 
discs. 



4.1 Case 1 

This case has parameters ^ = 5.58 {A = 4.0), ft = 0.1 and 
q = 0.005. We start measurements after the star is suffi- 
ciently resolved at the finest grid level. The radius (r,) is 
typically resolved by at least 20 cells. 



4-1.1 Self-limited spin up 

We first describe the early evolution 4 ^ t ^ 11, when the 
disc mass is insignificant compared to the stellar mass. Fig.[T] 
and Fig. [2] show snapshots of the density field and evolution 
of stellar properties, respectively. 

At t < 6.2, the star spins up with a rapid increase in 
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Figure 2. Evolution of stellar properties for Case 1 during the 
initial collapse phase. The top panel shows the spin to break-up 
frequency ratio (solid) and the ratio of kinetic to potential energy 
(dashed) . The bottom panel shows the orbital angular momentum 
(solid) and spin angular momentum (dashed). Note the different 
scales used on the left and right vertical axes. 



js, rJs/Slb and T/\W\ as it accretes material and angular mo- 
mentum. Because Q,s ~ js/ Sl, spin-up can occur because of 
decreasing 5**, but our EOS inhibits further collapse when 

p*. Hence, the spin-up is due to accretion, as expected 
for initial collapse. As it spins up, the star deforms into a 
bar-like object (Fig.[T] t — 6.66). The bar-shaped star exerts 
a positive gravitational torque on the surrounding material 
because it spins faster than the surrounding material, pro- 
ducing spiral arms in the latter (Fig.[TJ t = 7.23). This coun- 
teracts the increase in stellar spin angular momentum due 
to the accretion of material, resulting in an approximately 
constant js during 6.2 < i < 7.7. 

From t — 8 to t — 10 spiral arms are always present 
in the external fluid (e.g. Fig. [1] t = 9.25) and Fig. [2]shows 
that js is decreasing during this time interval. This indicates 
spin angular momentum loss due to negative torque exerted 
on the star by the prominent spiral arms. Thus the initial 
spin-up is limited by the increasing spin-down torque from 
the external fluid as the star becomes deformed into a non- 
axisymmetric, m — 2 object. Notice in Fig.[2]that jjoi ^ jjsl 
so the star is essentially fixed in the x — y plane in the 
inertial frame. However, note that the increase in jo at t = 
9.4 coincides with when js and fis/f^b stop decreasing. This 
suggests orbital motion hinders spin-down. 

Fig. E] compares the angular momentum fluxes at t = 
9.25, when js is decreasingj. The rise of gravity flux to- 
wards the average stellar surface, to aa ~ 0.5 at (r*) , 
implies some angular momentum loss due to gravitational 
torques. The figure shows that qa < around (r,) which 
indicates that large-scale advection is spinning up the star, 
but ac + Q.A ~ 0.2 > 0, implying that gravitational spin- 
down torques outweigh the effect of accretion, thereby pre- 
venting spin-up and leading to spin-down. The numerical 



^ Defining the external fluid relative to the star implies an in- 
direct potential associated with stellar motion, but this vanishes 
when performing azimuthal averages. 
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Figure 3. Case 1 early phase spin-down: non-dimensionalised, 
azimuthally averaged angular momentum fluxes due to the large- 
scale advection, self-gravity and Reynolds stresses. \R — Rt \ de- 
notes the cylindrical distance away from the star. The vertical 
line is (r») = 0.0146rc. This snapshot corresponds to the bottom 
right contour plot in Fig. [T] 
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Figure 4. Case 1 early phase spin-down. Azimuthally averaged, 
non-dimensionalised torque exerted by the star on its surround- 
ings. The region < O.OlSr'c corresponds to the star. 

viscosity is estimated to be ajv ~ 0.03 near p ~ p« and is 
therefore negligible. 

We also explicitly solve for the gravitational potential 
<1>, associated with p > p* and find a positive torque exerted 
by the star on the external fiuid, as shown in Fig. |4l This 
torque is concentrated near {r») . Together with the contour 
plots in Fig. [2] this is consistent with prominent m = 2 
spirals in the density field being responsible for spinning 
down the star. 



4-1.2 Long term evolution 

We ran Case 1 until the stellar spin approached a steady 
state. Its evolution is summarised in Fig. [5] Notice that 
the early phase discussed previously, when most variation 
in stellar spin occurs, coincides with the time when the pre- 
dicted disc is undefined. 

Between 11 < t < 26, fis/fib is highly variable but 
settles to ~ 0.5 at t = 34. Between 39 < t < 74 it decreases 
at a negligible rate compared to 7.7 J; t ^ 11- Most of the 
variation in r/| VF| also diminishes after t = 26. Notice that 
Q.s/0.h and T/\W\ become relatively constant as Ald/M, 
increases. 

The stellar spin angular momentum behaves similarly. 
There is no appreciable change in js relative to t < 11 after 
disc formation. However, between 11 ^ t < 20, the orbital 
angular momentum increases by two orders of magnitude 
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Figure 1. Case 1: contours of log (p/p, ) in tlie star's equatorial plane at four snapshots during the initial collapse. Thick lines delineate 
P = p.- 



(comparing jo in Fig. [5] and Fig. [S} from essentially zero to 
values comparable to js- at t ~ 20, jo — 0.056 and js — 
0.011. 

As Md/M, increases from ~ 0.3 to ~ 0.5, the orbital 
angular momentum jo undergoes large oscillations, which 
are comparable in magnitude to js. Indeed, we found that 
the star exhibits complex orbital motion, and its displace- 
ment from the box centre can be ~ 0.022rc, larger than its 
size ({r«) ~ 0.0146rc). Significant orbital motion coincides 
with disc formation. 

The evolutionary plots suggest that the disc inhibits 
changes to stellar spin. Instead, its interaction with the disc 
leads to orbital motion. This is consistent with a Fourier 
analysis of the surface density shown in Fig. [BlH- During 
the self-limited spin-up phase (t < 11), the m = 2 mode is 
dominant in the external medium to the star. The m — 1 



As the predicted disc does not develop until t ~ 20, strictly 
speaking amplitudes prior to this time are not disc modes but 



simply that of the external medium. 



mode over-takes m = 2 at f ~ 15, which is when jo begins 
to increase. After t ~ 20, the disc forms a strong m = 1 
asymmetry. The dominance of the m = 1 mode coincides 
with limited evolution in stellar spin, but large amplitudes 
in orbital motion. 

Early theoretical work shows that an m = 1 lopsided 
over-density causes st ellar orbital motion abo ut the sys- 
tem's centre of mass (| Adams et al.|[l98 9: Hcc mskerk et al.] 
Il992h . In those studies the star is treated as a point mass, 
so star-disc torques can only affect the star's orbital mo- 
tion. Although we model the star as a finite-sized poly- 
trope, and thus have an associated spin angular momentum, 
the angular momentum exchange between the star and the 
m = 1-dominated disc is still in orbital angular momentum. 
Changes to the spin angular momentum are apparently in- 
hibited. 

Fig. [7] shows the evolution of the dimensionless mass 
accretion rate. During the very early stages (t < 7), i.e. the 
initial spin-up phase, ^aim ~ 4 and roughly constant (recall 
the initialisation parameter ^ = 5.58). ^sim < C since the 
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Figure 5. Case 1: evolution of the disc-to-star mass ratio (top) 
and the same stellar properties shown in Fig. [2] (middle, bottom) 
over the entire simulation. Note the different scales used for the 
left and right vertical axes. 




Figure 6. Case 1: evolution of non-axisymmetric modes in the 
surface density of a region of cylindrical radius 0.25rc and thick- 
ness 0.025rc about the star. Fourier amplitudes were integrated 
from a cylindrical radius of 0.02rc to 0.25rc away from the star, 
then normalised by the axisymmetric amplitude. 



core has angular momentum so spherically symmetric col- 
lapse cannot occur. The sharp drop in ^sim in 7 < t < 10 
coincides with spin-down. As the star loses angular momen- 
tum to the immediate surroundings, the latter gains angular 
momentum, inhibiting accretion. However, once the star-on- 
disc torque is reduced (because the star has spun down and 
become more axisymmetric), material can again fall in and 
^sim rises {10 < t < 20). When the stellar spin is in steady 
state ^sini < but remains of order unity. The lower ^sim 
may result from the fact that at later times, the disc is lop- 
sided and is not well described by accretion from a spherical 
cloud onto an axisymmetric disc. 




Figure 7. Case 1: non-dimensionalised mass accretion rate. The 
large peak at t ~ 12 may be an artifact from numerical time 
derivatives. 



4-1.3 Structure and angular momentum transport at late 
stages 

In order to understand why spin evolution is inhibited on 
timescales beyond the early phase, we analyse the disc struc- 
ture while js and ^^/^h remain approximately constant. 
Fig. [8] and Fig. [9] respectively show the density field in real 
and Fourier space a.t t = 68.7. At the chosen snapshot, the 
star has size (r«) ~ 0.0116rc and the theoretical disc radius 
is Rk — O.lrc. Note that the lopsided disc is contained within 
Rk- We do not expect a sharp cut-off in Fourier amplitudes 
beyond R^ because the non-axisymmetric disc distorts the 
material beyond Rk. 

We see that m = 1 dominates the outer ~ 40% of the 
disc, corresponding to the lopsided density in the contour 
plot. The lopsided disc behaves like a binary companion to 
the star, causing the star to orbit about their common centre 
of mass. We do not expect, nor find, this motion to alter stel- 
lar spin. The fact that m = 1 is dominant towards the outer 
disc edge is favourable for inducing orbital motion because 
it acts like a lever arm. 

Close to the star, m = 2 dominates, but the contour 
plot shows spiral arms that are thin and much less promi- 
nent than those identified during initial collapse (Fig. [1]). 
The two arms also have different densities (reducing the am- 
plitude of m = 2 symmetry). Well beyond the disc edge Rk, 
m = 2 again dominates, but this low density core material 
here is not expected to have significant impact on the stel- 
lar spin because the total mass in this region is small, and 
because the star appears more axisymmetric as the m — 2 
component of its potential decays with distance. 

The limited spin evolution is due to ineffective gravi- 
tational coupling between the star and the disc. For a disc 
with m fold symmetry, a qualitative representation of its 
instantaneous surface density is 

E ~ ^ sin m(j!) -I- B cos m(/), (17) 

and for a star whose potential is dominated by k fold sym- 
metry, we can write 

^* A* sin k(j) + B* cos k(t>, (18) 

where A, B, A*, B, are real functions of radius. The instan- 
taneous torque per unit area exerted on the disc by the star 
is T = — Hence, for a disc dominated by the m = 1 
mode and a stellar potential with k = 2 symmetry (Fig. |8]), 
we expect {T}^ ~ 0. 
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Figure 10. Case 1: non-dimensionalised, azimuthally averaged 
angular momentum flux towards the end of the simulation. The 
vertical line is (r») . The advcctivc flux is ~ —0.5 and beyond the 
scale of this plot. 
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Figure 8. Case 1 towards the end of the simulation. The colour- 
bar indicates logp/p* in the star's equatorial plane. The inner 
black line indicates p = p* and the outer white circle indicates 
the disc edge ijj.. 
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Figure 9. Case 1 towards the end of the simulation; Fourier 
amplitudes of non-axisymmetric modes in surface density, corre- 
sponding to the snapshot in Fig.|8] The vertical line indicates the 
predicted disc radius R^. 

In other words, there is very httle torque if the two com- 
ponents do not share the same symmetry (m ^ fc). Phys- 
ically this is because there is no resonance, so for half the 
azimuth the torque is opposite sign to the other half, and 
averages to zero. We expect a much reduced self-gravity an- 
gular momentum flux across the steUar surface in such sit- 
uations. 

Fig. [To] shows azimuthally averaged angular momen- 
tum fluxes for the chosen snapshot. The self-gravity flux is 
ac — 0.05 near (r,) , and is an order of magnitude smaller 
than at i = 9.25 (Fig. [3]), indicating ineffective gravitational 
spin-down. This is consistent with the limited spin evolution 
observed during this time, and the fact that the star is much 
more axisymmetric than its bar-like shape at t < 11, thus 
more difficult to spin down. 

We find that limited spin evolution is also correlated to 



signiflcant motion of the star, induced by the m = 1 mode. 
We now examine a system with a smaller disc-to-star mass 
ratio, and consequently less orbital motion of the star. 

4.2 Case 2 

This case has m = 2 as the dominant non-axisymmetric 
mode throughout the simulation. Case 2 has parameters 
C = 2.74 {A = 2.8), h = 0.1 and q = 0.01. A lower ^ cor- 
responds to a smaller core and disc mass compared to Case 
1. Increasing q corresponds to a larger star initially. Again, 
we wait until the star is well-resolved at the finest grid level 
before making measurements. The star size (r*) is typically 
resolved by at least 30 cells. 

4-2.1 Star-disc evolution 

Fig. [11] shows the evolution of disc mass and stellar prop- 
erties. The phase f < 23 is similar to i < 11 of Case 1. 
There is an initial rapid increase in fis/nb, up to Q,b — 
0.5inb followed by rapid spin-down, abruptly stopping at 
Q.B — 0.498f2b. However, unlike Case 1, for t > 22> there is 
a noticeable, almost monotonic, decrease in f2s/f2b- Notice 
the plateau in spin frequency around t ~ 65, which coin- 
cides with increase in orbital motion (see below) when the 
disc becomes somewhat massive (M^ ~ O.OSM*). Like Case 
1 though, r/|W| remains fairly constant. The final spin fre- 
quency fis ~ 0.482f2b is not much lower than the maximum 
value attained during initial collapse, but the decreasing js 
suggests it may further spin-down if the simulation were 
continued. 

The evolution of the spin angular momentum is also 
much smoother compared to Case 1. Unlike the previous 
case, for t > 23, js decreases monotonically, reaching a value 
~ 10% lower than the maximum at t ~ 23. Here |jo| remains 
at least an order of magnitude smaller than js, with |jo| only 
growing when the disc develops (the latter was also seen for 
Case 1). 

At late times, the difference in stellar spin evolution 
between the present case and Case 1 is likely due to very 
different disc properties. In Case 2, Md/M, < 0.12, smaller 
than Case 1 by a factor of ~ 4. A disc with large M^^/Md is 
expecte d to be stable against t he growth of m = 1 pertur- 
bations (jHeemskerk et alii 19921 ). because there is insufficient 
disc mass to move the star. 

Fig. [12] shows the evolution of non-axisymmetric modes 
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Figure 13. Case 2: the dimensionless mass accretion rate. Note 
that our measurements begin after the central object is well re- 
solved. The initial phase observed in Case 1, where ^^i^ is roughly 
constant, is not present because for this simulation, the central 
object is not yet well resolved during these very early times. 



Yuan fc Cassenlll985h . This effect is absent in lKratter et alj 



2010l) . where ^ remains close to its initialisation value, be- 



cause in those simulations the central object was represented 
by a sink particle that was not capable of torquing the disc. 



Figure 11. Case 2: evolution of the disc-to-star mass ratio (top), 
the stellar spin and kinetic-to-potential energy ratio (middle) and 
the angular momenta (bottom) . Note the different scales used on 
the left and right vertical axes. 




Figure 12. Evolution of non-axisymmetric modes in the surface 
density of material external to the star. Fourier amplitudes were 
integrated over a cylindrical distance [0.04, 0.25]rc away from the 
star, then normalised by the axisymmetric amplitude. 



in surface density and is very different to Case 1. Here, the 
m — 2 mode remains dominant throughout the simulation. 
This is consistent with the idea that the m = 2 mode pro- 
vides the spin-down torque while the m = 1 mode produces 
orbital motion. 

Fig. [13] shows the measured mass accretion rate, ^aim 
is smaller than Case 1 because of the lower mass core. As 
before, there is a drop in accretion associated with self- 
limited spin. At late times, ^sim is roughly constant but is 
much smaller than ^. This could be because we have a disc 
which is receiving angular momentum from the spin-down 
of the central star. This can be thought of as an accretion 
disc with a positive torque applied at the inner boundary 



4-2.2 Disc structure and angular momentum flux 

We examine Case 2 at t — 59.7, during its long-term spin 
down. Fig. 1 141 and Fig. [15] show the density field in real and 
Fourier space, respectively. The average radius of the star is 
(r,) ~ 0.036rc and the theoretical disc edge is Rk — 0.106rc. 
The contour plot shows that the bulk of the disc is within 
Rk of the star, but it distorts the region beyond it. 

The contour plot clearly shows the star-disc system has 
m = 2 symmetry. Note that the relative Fourier amplitudes 
are quite different from Case 1, as the m — 1 mode is never 
dominant. Instead the m = 2 mode dominates most of the 
disc, and even the region beyond Rk- 

This explains why Case 2 experiences spin-down 
whereas Case 1 does not. In both cases the star has m — 2 
symmetry, but only in Case 2 does the disc also have a signif- 
icant m — 2 amplitude. Furthermore, the m — 1 amplitude 
in Case 2 is smaller than that of Case 1, consistent with 
limited orbital motion. The lack of a strong m — 1 mode 
appears favourable for spin-down via gravitational torques. 

Fig. [16] shows angular momentum fiuxes along the long 
and short axis of the star in its equatorial plane. In both 
cases, Reynolds stresses are negligible compared to gravity 
or advection. The gravity flux is large, with ac = 0{1). 
Note the sign of the fluxes depends on the azimuthal direc- 
tion. Averaging the fiuxes over the non-axisymmetric radius 
where p = p*, we found ao ~ 0.3, an = 0(10"^) and the 
large-scale advective flux aA ~ —0.1. The numerical vis- 
cosity is Qjv = O(10~^). The total non-dimensional fiux is 
~ 0.2, consistent with the observation of spin angular mo- 
mentum loss. 

We computed the potential $, to find the torque ex- 
erted by the star on the surrounding disc, shown in Fig. 1171 
The torque becomes negative and large in amplitude towards 
0.04rc, or approximately the stellar surface. This large neg- 
ative torque might lead to the impression that the disc is 
spinning up the star, contradicting the observed spin-down, 



10 Lin, Krumholz & Kratter 

-4.0 -3.2 -2.3 -1.5 -0.7 0.7 ' .0 
0.1 6 I 




0.16 



Figure 14. Case 2: snapshot of log (p/pt) in the star's equatorial 
plane. The inner black line indicates p = pt, and the outer white 
line indicates iJ^.. 
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Figure 15. Case 2: Fourier amplitudes of non-axisymmetric 
modes in the surface density corresponding to Fig. 1141 The pre- 
dicted disc edge is marked by the vertical line (right). The star 
has an average radius (r*) ~ 0.036rc. 



but care must be taken when interpreting azimuthally av- 
eraged, one-dimensional profiles such as Fig. [iTl Since the 
star-disc interface is non-axisymmetric, there is no single 
value of radius to represent the star-disc interface. However, 
the azimuthally-averaged torque is positive around 0.06rc 
from the star, which is certainly in the disc region (see Fig. 
[14] where this region has typical densities much lower than 
the transition density). This implies a net loss of spin angu- 
lar momentum from the star to this disc region. 



5 DISCUSSION 

In our simulations, we observe an upper limit to stellar spin 
that is below break-up speed. This limitation occurs during 
the initial spin-up of the star. This relies on the ability of 
the central object to deform into a non-axisymmetric shape 
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Figure 16. Case 2: non-dimensionalised angular momentum 
fluxes during the slow spin-down phase. Top: fluxes along the 
<j) = it/4 azimuth (approximately the long axis of the star. Bot- 
tom: fluxes along the <j) = 3n/4 (approximately the short axis of 
the star). In each, the vertical line indicates the surface of the 
star along that direction. Note the different scales used for the 
two panels. 
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Figure 17. Case 2: azimuthally averaged torque per unit area 
exerted on the disc by the star, during the slow spin-down phase. 



and lose angular momentum to its surroundings by exerting 
gravitational torques on the external medium. 

Consider a rotating, axisymmetric self-gravitating fluid 
body of mass AI, angular momentum J and radius a in the 
plane perpendicular to the spin axis. We define the dimen- 
sionless angular momentum as 

(19) 



7^ 



Now, specialise to a Maclaurin spheroid with angular mo- 
mentum J = (2/5)a^cjM, where ui is the uniform spin an- 
gular frequency. The dimensionless parameter TZ for the 
Maclaurin spheroid is 



GM 

r,3 



7^^ 



(20) 



If TZ-Mac ^ .278, triaxial equilibri um with uniform rotation 
is possible (|Yuan fc CassenI Il985l). The quantity in square 
brackets is equivalent to our spin parameter Qg/Qh- If we 
increase 7?.Mac from a small value, the Maclaurin spheroid 
will exceed the critical value above before it reaches break-up 
speed (at which point 7?.Mac = 0.4). That is, it may become 
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non-axisymmetric before flying apart. Obviously, our cen- 
tral object is not a Maclaurin spheroid, but inserting the 
observed maximum spin for Case 1, f^s/ilb ~ 0.55, yields 
7?.Mac = 0.22, which is not too different from 0.278. 

This analysis makes an important point which, together 
with our our simulations, suggests that the maximum stel- 
lar spin may be determined largely by the critical value of 
TZ (or Qs/^h), beyond which instabilities in the rotating 
object produce non-axisymmetric deformations. As soon as 
this occurs, the object may torque the surrounding disc and 
lose angular momentum to it. The transition from Mac lau- 
rin spheroids to Jacobi ellipsoids (|Chandrasekhail|l969i ) at 
TZ = 0.278 is one example, but in general the critical value 
will depend on the pr operties of the central object, such 
as its equation of state (|Ostriker fc Bodenheimenl 19731) . We 
may think of objects with lower critical values of TZ as being 
easier to deform, and will likely be limited to smaller max- 
imum spin rates. For example, our compressible star with 
smaller 7?.Mac is easier to deform than the incompressible 
Maclaurin spheroid. 

In our simulations we treat the star as a simple n — 3/2 
polytrope, but this is obviously an oversimplification of true 
protostellar structure, which will in general depend on its 
mass and evolutionary state. Even if we were to include a 
more realistic treatment of the EOS for stellar material, our 
limited resolution would make it impossible to study the 
star's deformation in anything but a very schematic man- 
ner. Such a study is beyond the scope of this paper, but due 
to this limitation we should be cautious about putting too 
much weight on the exact numerical values we compute for 
the gravitational spin-down torque. Our general conclusions 
that the initial spin-up is self-limited to well below breakup, 
and that a prominent m = 2 mode in the surrounding disc is 
capable of providing a further, long-term spin-down torque, 
should be robust. In particular, it is important to note that 
the argument we have given above is completely dimension- 
less, and does not depend in any way on the true physical 
sizes of protostars or their discs. 

A related issue is the degree of deformation needed 
for effective angular momentum exchange between the star 
and its surroundings. In our simulation, the density changes 
smoothly across the star-disc interface, whereas in reality ra- 
diation from both the accretion shock and the stellar surface 
reduce the entropy of stellar material well below that of disc 
material. This causes stellar matter to be far denser than 
matter in the disc immediately outside the star, an effect 
that we miss because our simulations are non-radiative. A 
larger density contrast would make the star's gravitational 
field stronger, and would therefore increase the effective- 
ness of gravitational angular momentum transport relative 
to that in our simulations. As a result, real stars proba- 
bly spin down more effectively than we have found, and re- 
quire less non-axisym metry to achieve the same torque (c.f. 
I Yuan fc Cas"senlll985h . Thus a more realistic treatment of 
stellar structure may actually strengthen the effect we have 
identified. 

5.1 Numerical issues 

The issue of numerical viscosity was not discussed in detail, 
which could be important for numerical studies of rotating 
fiuids in a Cartesian box. Angular momentum conservation 



can be violated near the centre of the rotating fiuid. How- 
ever, we have re-calculated spin angular momenta excluding 
the innermost few cells and found negligible effect on their 
evolution. In addition, lower resolution runs were performed 
which yield similar angular momenta evolution. We are thus 
confident that the observed correlation between spin evolu- 
tion and evolution of disc modes are physical results. 

We have focused on only two simulation runs. Although 
insufficient to predict stellar spin evolution as a function of 
the problem parameters, our results provide an useful guide 
to future numerical studies. We expect discs that develop 
prominent global m — 2 spirals to be most effective at spin- 
ning down the star through gravitational torques, because 
the star will usually have the same symmetry. 

5.2 Additional caveats 

It should be noted that the physics in our problem setup is 
highly simplified. The dynamical range covered in our sim- 
ulations does not correspond to an entire, realistic star-disc 
system. This discrepancy could potentially limit the appli- 
cability of our results to realistic star-disc systems. 

One issue is that if the star-on-disc torque is most sig- 
nificant close to the stellar surface, then such torques will 
be ineffective if there exists a significant gap between the 
star and the inner disc edge. A gap may be caused by mag- 
netic fields for low ma ss stars. However, using the model of 
iMatt fc Pudrit3 (|2005l ). one finds that for massive stars the 
truncation radius is actually inside the star, implying no gap 
( Rosen et al. 2011,) . Gravitational torques may then play a 
role. In addition, if gravitational torques are a result of high 
protostar spin rate, which enables triaxiality, then material 
must have fallen onto the star, adding to its angular mo- 
mentum in the first place. The smooth star-disc boundary 
in our simulations is then self-consistent with limiting spin 
via gravitational torques in the early stages of collapse. 

Finally, although we have demonstrated the existence 
of star-disc torques, we caution that the exact magnitudes 
of the torques will likely depend on the structure of the disc, 
which can in turn be modified by magnetic fields and radi- 
ation. We have not included these effects, and in their pres- 
ence the true spins of stars undergoing gravitational torque- 
regulated spin-up as they form will likely be quantitatively 
different than what we have found, but our general conclu- 
sions should hold. 



6 CONCLUSIONS 

We performed hydrodynamic simulations of the collapse of 
an isothermal cloud leading to the formation of a disc with 
a dense central object. We focused on the evolution of the 
central object's spin and its relation to disc properties. 

We find that the initial spin-up of the central object 
does not exceed ~ 50% of its break-up speed because in- 
creasing spin also increases its deformation into a bar-like 
object (possibly through an instability), on which the exter- 
nal material exerts a negative torque. Spin-down occurs over 
a period that is short compared to star formation timescales. 

We also find that when the surrounding disc develops a 
dominant m = 1 disturbance, the central object's spin evo- 
lution is inhibited, and large orbital motions are induced 
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by the disc. Our experiments suggest, perhaps counter- 
intuitively, that massive discs are less able to provide long 
term gravitational spin-down torques on a star because such 
discs are prone to developing m = 1 lopsided modes. 

Although spin-down was observed over the simulation 
timescale for systems with m = 2 symmetry, the spin-down 
rate is small compared to that during the earliest phase 
of collapse. The precise value of this upper limit may de- 
pend on the internal structure of the central object in ways 
that we have not explored. In general we expect that struc- 
tures that are easier to deform, in the sense that a smaller 
spin-to-break-up frequency ratio is needed to allow non- 
axisymmetry , will yield smaller maximum spin rates than 
structures that require larger spin-to-break-up frequencies 
to become non-axisymmetric. 

An upper limit to spin rate has important implications 
for stellar evolution. Strong rotational mixing may occur 
for large spi n rates and allow the star to bypass the re d 
giant phase (|Wooslev fc Hegeill2006l: lYoon fc Langeij|2005h . 



The critical spin rate found by I Wooslev fc Hegeii is :^_^0% 



Ekstrom et all (|2008l ) 



of break-up. As another example 
found that high rotation (~ 40% — 70% break-up) can 
increase metal production in initially metal-free stars. 
However, if gravitational torques limit stellar spin to ~ 50% 
of break-up, strong rotational mixing or increased metal 
production may be difficult to achieve because once the star 
reaches the main sequence, stellar winds provide further 
spin-down. A maximum stellar spin set by gravitational 
torques can be used to constrain the parameter space for 
stellar evolution calculations. 
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